home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / vfft / inrlvfft.cxx < prev    next >
Encoding:
Text File  |  1992-06-01  |  4.6 KB  |  151 lines
  1. /* inrlvfft.cxx
  2. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  3.  
  4. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  5. Permission is granted to reproduce this software for non-commercial
  6. purposes provided that this notice is left intact.
  7.  
  8. It is acknowledged that the U.S. Government has rights to this software
  9. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  10. and the University of California.
  11.  
  12. This software is provided as a professional and academic contribution
  13. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  14. with no warranties of any kind whatsoever, no support, no promise of
  15. updates, or printed documentation. By using this software, you
  16. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  17. University of California shall have no liability with respect to the
  18. infringement of other copyrights by any part of this software.
  19.  
  20. For further information about this notice, contact William Johnston,
  21. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  22. (wejohnston@lbl.gov)
  23.  
  24. For further information about this software, contact:
  25.     Jin Guojun
  26.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  27.     g_jin@lbl.gov
  28.  
  29. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  30.     */
  31. dimen2size = dimen1len * row << 1;
  32.  
  33. if (uimg.in_form==IFMT_DVFFT2D || uimg.in_form==IFMT_DVFFT3D)    {
  34. double    *otemp,
  35.     *obuf = nzalloc(fsize, sizeof(*obuf), "fobuf"),
  36.     *ibuf = nzalloc(dimen2size*(dimens?1:frm), sizeof(*ibuf), "fibuf");
  37.  
  38.     if (uimg.o_form==IFMT_BYTE)
  39.         otemp = nzalloc(fsize, 1, "otemp");
  40.     else    otemp = obuf;
  41.  
  42.     DBwork_space_init(MAX(MAX(cln, row), frm));
  43.     if (!dimens && frm>1)    {
  44.     int    r;
  45.     register DBCOMPLEX    *cmptr, *cp;
  46.     register int    c, n=frm, dimen1size = dimen1len;
  47.  
  48.         fsize = n * dimen2size;
  49.         i = fread(ibuf, sizeof(*ibuf), fsize, in_fp);
  50.         if (i != fsize)
  51.             syserr("3d[%ld] read %d", fsize, i);
  52.  
  53.         load_DBw(frm);
  54.  
  55.         for (r=0; r<row; r++)
  56.             for (c=0; c<dimen1size; c++) {
  57.             cp = (DBCOMPLEX*)ibuf + r*dimen1size + c;
  58.             cmptr = DBcmpx_in;
  59.             for (i = 0; i < n; i++, cp+=dimen2size>>1) {
  60.                 c_re (cmptr [i]) = c_re(*cp);
  61.                 c_im (cmptr [i]) = -c_im(*cp);
  62.             }
  63.  
  64.             DBFourier(cmptr, n, DBcmpx_out);
  65.  
  66.             cp = (DBCOMPLEX*)ibuf + r*dimen1size + c;/* store back to src */
  67.             cmptr = DBcmpx_out;
  68.             for (i = 0; i < n; i++, cp+=dimen2size>>1)
  69.                 *cp = cmptr[i];
  70.             }
  71.     }
  72.  
  73.     fsize = row * cln;
  74.     for (f=0; f<frm; f++)    {
  75.         if (dimens)    {
  76.         register DBCOMPLEX*    cmptr = (DBCOMPLEX*) ibuf;
  77.         i = fread(ibuf, sizeof(*ibuf), dimen2size, in_fp);
  78.         if (i != dimen2size)    syserr("2d[%d] read %d", dimen2size,i);
  79.         i >>= 1;
  80.         while (i--)
  81.             c_im(cmptr[i]) = -c_im(cmptr[i]);
  82.         }
  83.  
  84.         DBvrft2d(ibuf + (dimens ? 0 : f*dimen2size), cln, row, obuf);
  85.         if (uimg.o_form != IFMT_FLOAT)
  86.         dtob(obuf, otemp, fsize, 1./(fsize*(uimg.in_form&1 ? frm:1)));
  87.         else
  88.         dtof(obuf, otemp, fsize, 1./(fsize*(uimg.in_form&1 ? frm:1)));
  89.         i = fwrite(otemp, uimg.pxl_out, fsize, out_fp);
  90.         if (i != fsize)    syserr("3d[%d] write %d", fsize, i);
  91.     }
  92. }
  93. else    {
  94. float    *otemp,
  95.     *obuf = nzalloc(fsize, sizeof(*obuf), "fobuf"),
  96.     *ibuf = nzalloc(dimen2size*(dimens ? 1 : frm), sizeof(*ibuf), "fibuf");
  97.  
  98.     if (uimg.o_form == IFMT_BYTE)
  99.         otemp = nzalloc(fsize, 1, "otemp");
  100.     else    otemp = obuf;
  101.  
  102.     work_space_init(MAX(MAX(cln, row), frm));
  103.     if (!dimens && frm>1)    {
  104.     int    r;
  105.     register COMPLEX    *cmptr, *cp;
  106.     register int    c, n=frm, dimen1size = dimen1len;
  107.  
  108.         fsize = n * dimen2size;
  109.         i = fread(ibuf, sizeof(*ibuf), fsize, in_fp);
  110.         if (i != fsize)
  111.             syserr("3d[%ld] read %d", fsize, i);
  112.  
  113.         load_w(frm);
  114.  
  115.         for (r=0; r<row; r++)
  116.             for (c=0; c<dimen1size; c++) {
  117.             cp = (COMPLEX*)ibuf + r*dimen1size + c;
  118.             cmptr = cmpx_in;
  119.             for (i = 0; i < n; i++, cp+=dimen2size>>1) {
  120.                 c_re (cmptr [i]) = c_re(*cp);
  121.                 c_im (cmptr [i]) = -c_im(*cp);
  122.             }
  123.  
  124.             Fourier(cmptr, n, cmpx_out);
  125.  
  126.             cp = (COMPLEX*)ibuf + r*dimen1size + c;/* store back to src */
  127.             cmptr = cmpx_out;
  128.             for (i = 0; i < n; i++, cp+=dimen2size>>1)
  129.                 *cp = cmptr[i];
  130.             }
  131.     }
  132.  
  133.     fsize = row * cln;
  134.     for (f=0; f<frm; f++)    {
  135.         if (dimens)    {
  136.         register COMPLEX*    cmptr = (COMPLEX*) ibuf;
  137.         i = fread(ibuf, sizeof(*ibuf), dimen2size, in_fp);
  138.         if (i != dimen2size)    syserr("2d[%d] read %d", dimen2size,i);
  139.         i >>= 1;    /* size float -> complex    */
  140.         while (i--)
  141.             c_im(cmptr[i]) = -c_im(cmptr[i]);
  142.         }
  143.  
  144.         vrft2d(ibuf + (dimens ? 0 : f*dimen2size), cln, row, obuf);
  145.         if (uimg.o_form != IFMT_FLOAT)
  146.         ftob(obuf, otemp, fsize, 1./(fsize*(uimg.in_form&1 ? frm:1)));
  147.         i = fwrite(otemp, uimg.pxl_out, fsize, out_fp);
  148.         if (i != fsize)    syserr("3d[%d] write %d", fsize, i);
  149.     }
  150. }
  151.